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Abstract. 

We review the classification of all the 36 possible gap-opening instabilities in 
graphene, i.e., the 36 relativistic masses of the two-dimensional Dirac Hamiltonian 
when the spin, valley, and superconducting channels are included. We then show 
that in graphene it is possible to realize an odd number of Majorana fermions 
attached to vortices in superconducting order parameters if a proper hierarchy of 
mass scales is in place. 
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1. Introduction 

The discovery [T] that it is possible to peel and isolate individual atomic layers of 
graphite, i.e., graphene, has led to an explosion of experimental and theoretical 
works and ideas from exotic physics to real material applications [2J. Graphene is 
a material with remarkable physical properties, many of which are consequence of its 
band structure at charge neutrality, which is characterized by two Fermi points in the 
Brillouin zone. An excitation around those Fermi points, being linearly proportional 
to its momentum, resembles the dispersion relation of massless relativistic particles. 
The low energy theory is then well described by a relativistic Dirac Hamiltonian in 
two dimensions [3]. Massless Dirac fermions describe graphene, but one may wonder 
what kinds of mass gaps could be induced. This question is rather important for device 
applications, as the presence of a gap would make graphene behave as a semiconductor, 
just like silicon. But the question is also interesting for fundamental physics reasons, 
as we discuss below. 

A mass in the Dirac equation can be viewed as an order parameter, a bosonic 
field. If the mass is generated by spontaneous breaking of a symmetry - the Higgs 
mechanism in graphene - there can be spatio-temporal fluctuations in the mass order 
parameter. In particular, there could be topological defects in the order parameter: 
domain walls, vortices, or hedgehogs, for instance. Alternatively, the masses could 
be induced externally, for example if the gaps are attained by placing graphene on 
a certain substrate, such as one where there is a difference in the potential seen by 
the two atoms in the unit cell of graphene, or if the substrate is a superconductor. 
Topological defects could also be present in this case where the gaps are externally 
induced. 

Topological defects in an order parameter can lead to zero modes in the fcrmionic 
spectrum, i.e., zero energy solutions lying in the middle of the gap [3]. The physical 
consequences of these zero modes are rather remarkable. When charge is conserved, 
the zero modes bind a fraction of an electron charge [5], while in superconducting 
systems, they bind charge neutral Majorana fermions [6]. 

In this paper, we classify all possible 36 competing orders of a Dirac Hamiltonian 
represented by 16-dimensional Dirac matrices that encode the spin, valley, and 
superconducting channels. We also discuss the simpler cases where only spin and valley 
degrees of freedom (no superconductivity), or valley alone (spinless electrons) are 
considered. These simpler cases serve as a warm up exercise in gaining familiarity with 
the classification construction as the build up of the increasingly larger representations 
is carried out. 

We then show how this classification can be applied to the problem of selecting 
an odd number of Majorana fermions to bind to superconducting vortices. That 
the number of zero modes must be odd is important for applications to topological 
quantum computing, that we shall discuss in more detail below. It is possible to achieve 
an odd number of Majorana fermions using surface states of topological insulators [?]■ 
But, naively, Majorana zero modes bound to superconducting vortices in graphene, 
first found by Ghaemi and Wilczek [6] in this context, would come in quadruplets 
because of the valley and spin degrees of freedom. We show otherwise. If there is a 
proper hierarchy of mass scales, one can tune selectively the number of Majorana zero 
modes in graphene from 4 — >■ 3 — > 2 1 — > 0. 

We aim in this paper at a pedagogical description of the mass classification 
scheme, and a brief discussion of how the number of zero modes can be tuned by 
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changing the relative strength of masses. The original classification of masses in 
graphene was presented in [5], while the tuning of the number of Majorana modes 
was carried out in [5]. We point the reader to those papers for details beyond this 
review. 



2. Classification of masses in graphene 

Let us first of all define what we mean by a mass in graphene. Say we take generically 
a two-dimensional Dirac Hamiltonian 

H = p x a 1 + p y a 2 + m M (1) 

where p x y stand for momentum operators in the x, y-directions, and a 1 2 are some 
(generically) D-dimensional matrices satisfying {a^a^} = 2 6^. We would like to say 
that m is a mass scale, and that M is a mass matrix. For that to be the case, the 
matrix M has to satisfy certain commutation relations with the matrices a 1 2 . The 
relations are seen if we square the Dirac-type Hamiltonian: 

H 2 =p 2 x +pl + m 2 M 2 +mp x {a 1 ,M} + mp y {a 2 ,M}, (2) 

which if 

M 2 = 1 and {a 12 ,M} = 0, (3) 

yields 

H 2 =vl +Py +m 2 . (4) 

Thus a mass matrix is one that satisfies the relations in Eq. ([3]) and leads to the 
dispersion for a massive Dirac particle e p = ±\/IpI 2 + m2 ■ 

Let us then count the number of such matrices for cases of increasing complexity, 
starting with the case of spinless fermions, marching along to the most general, 
when we take into account both spin and particle-hole gradings needed to discuss 
superconducting graphene. 



2.1. Spinless case - A- dimensional Dirac matrices 

Let us write the 4x4 Dirac matrices that describe spinless electrons in graphene in 
the Weyl representation: 

= ^3 ® r f , (3 = ( ° Q q° J = o x ® Tq , (5) 

where the 2x2 unit matrix r and the three Pauli matrices r x , r 2 , and t 3 act on 
the sublattices indices (A, B), while the 2x2 unit matrix er and the three Pauli 
matrices a 1 , cr 2 , and cr 3 act on the valley indices (+,—). In this representation, the 

four-component spinor is given by \E , t = (^<p^ +B , V'+A' ^— A' V'-s) • 

We take the matrices a 1 2 to construct the kinetic part of the Hamiltonian, such 
that the gapless case can be written as 

=Px a i +P y "2- (6) 
Now, the possible masses in the 4-dimensional representation of the Dirac Hamiltonian 
corresponds to all matrices M of the form X^ v = cr^® t v (other than a 1 = X 31 and 
a 2 = X 32 ) such that {X^^a^} = and {X^,a 2 } = 0. The list of such matrices has 
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4 elements: X 33l X w , X 20 , and A" 03 . These are then the 4 possible mass terms in 
spinlcss graphene. 

Physically, we can identify these four mass matrices as follows. One perturbation 
that opens a gap 2/i s is a staggered chemical potential, taking values +fj, s and — fi B 
in the two sublattices A and B of graphene. This is the mass term added by 
Semenoff [10], and it corresponds to A 33 . A second mass gap 2 1 77 1 arises by adding 
directed next-nearest-ncighbor hopping amplitudes in the presence of fluxes, but such 
that no net magnetic flux threads a hexagonal Wigncr-Scitz unit cell of graphene. 
This perturbation, introduced by Haldane breaks time- reversal symmetry (TRS) 
and corresponds to X 03 . Finally, a real- valued modulation of the nearest-neighbor 
hopping amplitude with a wave vector connecting the two Dirac points (i.e., a Kekule 
dimerization pattern for graphene) also opens a gap 2|A| [S]. This real- valued 
modulation of the nearest-neighbor hoppings is parametrized by the complex order 
parameter A = Re A + i Im A whose phase controls the angles of the dimerization 
pattern. This complex order parameter translates into two real masses Re A and 
ImA, corresponding respectively to X w and X 20 , bringing the total number of real- 
valued masses to four. 

We can identify the microscopic origin of these 4 masses according to the 
symmetries that they respect (or break) at the level of the Dirac equation. To this 
end, let us look into two symmetries that the Dirac Hamiltonian H may or may not 
possess: 

• Time-reversal symmetry (TRS) is satisfied when 

X 11 H*(-p)X n = H(p). (7) 

• Sub-lattice symmetry (SLS), also referred to as chiral symmetry in the literature, 
is satisfied when 

A 33 H{p) X 33 = -H{p) . (8) 

This symmetry operation corresponds to flipping the sign of the wavefunction in 
one of the sublattices but not on the other (here, we flip the sign on sublatticc 
A). 

The Kekule dimerization pattern corresponds to a spatial modulation of the 
hopping matrix elements between sublattices A and B, and therefore changing the sign 
of the wavefunction on one sublattice would reverse the overall sign of the Hamiltonian. 
And indeed, one verifies that masses of the form X 10 and X 20 satisfy Eq. ((8]). The 
hopping amplitudes in the Kekule dimerization pattern are all real, and therefore they 
respect TRS, which can also be checked via Eq. ([7]). 

Both the staggered chemical potential or Semenoff mass X 33 and the Haldane 
mass X 03 break SLS, as they involve couplings between sites in the same sublattices. 
However, the staggered chemical potential respects TRS, while the Haldane mass 
breaks it, as can be checked using Eqs. (0 and (0). 

We shall introduce here a terminology that will be useful later on when we 
consider systems with larger size representations, once spin degrees of freedom and 
superconductivity are considered. For instance, the Kekule dimerization pattern 
corresponds to a valence-bond solid (VBS) order parameter by analogy to the 
terminology used for quantum dimer models. A VBS order picks up a microscopic 
orientation that translates into a complex-valued order parameter in the continuum 
limit. Hence, we shall distinguish between the real (Re VBS) and imaginary (ImVBS) 
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X w RcVBS True True 

X 20 ImVBS True True 

X 33 CDW True False 

X 03 QHE False False 



Table 1. The 4 mass matrices that can be added to the massless Dirac 
Hamiltonian K, from Eq. (J6J arc of the form = a eg) t v and anticommutc 

with /Cq . Each mass matrix can be assigned an order parameter for the underlying 
microscopic model. Each mass matrix preserves or breaks time-reversal symmetry 
(TRS), see Eq. Q and sublattice symmetry (SLS), see Eq. J5J. 



parts of the VBS, and here they correspond to the Re A and ImA of the Kckulc 
distortion. Now, any mass matrix that does not satisfy the SLS defined in Eq. ((SJ 
corresponds to a microscopic order parameter for which the fermion bilinear has the 
two lattice fermions sitting on the same sublattice. This is the case of the mass 
associated to the staggered chemical potential, which we may identify with a charge- 
density wave (CDW). Finally, the Haldane mass implies a quantum Hall effect (QHE). 
This nomenclature, along with the classification of the 4 masses according to TRS and 
SLS, are summarized in Table [T] 



2.2. Spinful case - 8- dimensional Dirac matrices 

Having warmed up with the 4x4 representations of the Dirac matrices for the simpler 
spinless case, we now construct the representations of all the masses for the case of 
graphene with the spin degrees of freedom included, but still no superconductivity 
considered. 

To represent the single particle Hamiltonian /C, utilize the 64 8-dimensional 
Hcrmitean matrices 

where pL l 2 3 = 0,1,2,3. Here, we have introduced the three families , cr^ , and 
of unit 2x2 and Pauli matrices that encode the spin- 1/2, valley, and sublattice 

degrees of freedom of graphene, respectively. 

The masses should be added to the massless Dirac Hamiltonian 

= Px a i +P V a 2, where a x = X 031 and a 2 = X Q32 . (10) 

The possible masses in the 8-dimensional representation of the Dirac Hamiltonian 
correspond to all matrices M of the form (other than a 1 and a 2 ) such that 

{ x ^^ 2 ^ a i\ = and { x MlM3 n 3 > a 2} = °- 0ne can carr y out the exercise of finding 
such matrices, obtaining the 16 matrices listed in Tabic 

We can classify these matrices according to the symmetries that they respect (or 
break). The symmetries that the Hamiltonian H may possess are as follows: 

• Time-reversal symmetry (TRS) is satisfied when 

X 211 H*(-p)X 211 =H(p) . (11) 
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• Sub-lattice symmetry (SLS) is satisfied when 

X 033 H(p) X 033 = -H(p) . (12) 

• Spin-rotation symmetry (SRS) is satisfied when 

[X 100 ,H(p)) = [X 200 , H(p)} = [X 300 , H(p)} = . (13) 

Notice that once we introduce spin degrees of freedom, we can decide whether the 
Hamiltonian is spin rotational invariant or not. Also, we can introduce terminology 
similar to the VBS and CDW cases we used in labeling the masses for the spinless case. 
Masses that do not satisfy the SLS defined in Eq. (fT2]l correspond to a microscopic 
order parameters for which the fermion bilinear has the two lattice fermions sitting 
on the same sublattice. Microscopic examples are, in addition to the CDW already 
previously encountered, the spin-density waves (SDW) such as Neel ordering, orbital 
currents leading to the quantum Hall effect (QHE), and spin-orbit couplings leading to 
the quantum spin Hall effect (QSHE). Whenever the instability can have a direction 
associated to it in internal spin space, we add the corresponding directional subscript 
x,y,z. 

The first 4 masses listed in Table [5] are physically the same as those 4 listed 
in Table [T] for the case of spinless electrons. These 4 masses correspond to order 
parameters in the charge sector. The next 12 masses correspond to some form of 
magnetic order. The simpler are the Neel^ y z order parameters along the three 
directions. The Neel states are the SDW order associated to fermion bilinears at the 
same lattice site. The RcVBS^ y z and ImVBS^, z are their counterparts where the 
fermion bilinears are defined on the bonds instead of the sites. Finally, the QSHE^. z 
correspond to the quantum spin Hall effect discussed by Kane and Mclc in Ref. [l"2j . 



2.3. Fully general case of single-layer graphene - 16- dimensional Dirac matrices 

To describe all symmetry-breaking instabilities with a local order parameter in 
graphene we consider the Bogoliubov-dc Gcnnes (BdG) Hamiltonian 

H BdG = \ [ d 2 r¥ Hi> (14) 



2 

where ^> is the 16-component Nambu spinor 

",t „?,t 



* = i>\, r v % ) (15) 

and ip s —^ ^ is a 4-component fermion annihilation operator that accounts for the 2 
valley and the 2 sublattice degrees of freedom. The kernel of the BdG Hamiltonian 
has the block structure 

H=(y u $ ) (i6) 

where the 8x8 blocks TL pp and H ph act on the combined space of valley, sublattice, 
and spin degrees of freedom, and represent the normal and anomalous part of the BdG 
Hamiltonian, respectively These blocks satisfy 

U PP ='Hp P (Hermiticity) , Hp h = — H ph (Fermi statistics). (17) 

To represent the single particle Hamiltonian H, define the 256 16-dimensional 
Hcrmitean matrices 

= P Ml ® ® <V 3 ® T A* 4 ( 18 ) 
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Mass matrix 


Order parameter 


TRS 


SRS 


SLS 


^010 


ReVBS 


True 


True 


True 


V 

^020 


1111 V -DO 


True 


True 


True 


^033 


CDW 


True 


True 


False 


^003 


QHE 


False 


True 


False 


^110 


ReVBS^ 


False 


False 


True 


^210 


ReVBS y 


False 


False 


True 


-^310 


ReVBS 2 


False 


False 


True 


^120 


ImVBS^ 


False 


False 


True 


^220 


ImVBS y 


False 


False 


True 


^320 


ImVBS^ 


False 


False 


True 


^"l03 


QSHE^ 


True 


False 


False 


^203 


QSHE y 


True 


False 


False 


^303 


QSRE Z 


True 


False 


False 


^133 


Neel x 


False 


False 


False 


^233 


Neel y 


False 


False 


False 


^333 


Neel z 


False 


False 


False 



Table 2. The 16 mass matrices that can be added to the massless Dirac 
Hamiltonian K. Q from Eq. H1UII are of the form X . = s ^ ® a v % t a and 
anticommute with K Q . Each mass matrix can be assigned an order parameter for 
the underlying microscopic model. The latin subindcx of the order parameter's 
name corresponds to the preferred quantization axis in SU(2) spin space. Each 
mass matrix preserves or breaks time-reversal symmetry (TRS), see Eq. Ijllll . 
spin- rotation symmetry (SRS), see Eq. 1121 1 . and sublatticc symmetry (SLS), see 
Eq. (13). 

where ^ 2 3 4 = 0, 1, 2, 3. Here, we have introduced the four families , , , and 
of unit 2x2 and Pauli matrices that encode the particle-hole (Nambu), spin- 1/2, 
valley, and sublattice degrees of freedom of graphene, respectively. 

The Dirac kinetic energy /C of graphene that accounts for the BdG block 
structure (I16|) is given by 

K, = p x ati +p y a 2 , where a x = X 0031 and a 2 = X 3032 . (19) 
There are 64 = 4 x 16 mass matrices (i.e., X,, ,, ,, ,, that anticommutes with 
/Cq). Of these 64 mass matrices, only 36 satisfy the condition 

-^1000 -^HyH^HsHi ^1000 = ~-^/j. 1 n 2 n 3 n i (20) 

for particle-hole symmetry (PHS) and are thus compatible with the symmetry 
condition (p x <E) s <E> a ® t ^)* = on the Nambu spinors [i.e., compatible with Eq. 
(1171) ]. All mass matrices with PHS are enumerated in Table [3] 

All 36 mass matrices from Table [3] can be classified in terms of the following 
(microscopic) 3 symmetry properties. 
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• A BdG Hamiltonian has time-reversal symmetry (TRS) when 

X 02n H*(-p)X 0211 = H(p). (21) 

• A BdG Hamiltonian has sublattice symmetry (SLS) when 

X 0Q33 H(p)X 0033 = -H(p). (22) 

• A BdG Hamiltonian has SU(2) spin rotation symmetry (SRS) when 

[*3ioo. H(jp)\ = [A 0200 , H(p)} = [X 3300 ,H(p)} = 0. (23) 

Identifying the physical meaning of the masses is done in a similar way as 
explained in the simpler cases discussed previously (the spinless and spinful cases 
without superconductivity). Below we present the rational for the terminology in full 
generality 

The microscopic order parameter corresponding to a mass matrix satisfying the 
SLS defined in Eq. (f2"2"j) is a non-vanishing expectation value for a fermion bilinear with 
the two lattice fermions residing on the ends of a bond connecting sites in opposite 
sublattices. We shall say that such a mass matrix is associated to a valence-bond solid 
(VBS) order parameter in analogy to the terminology used for quantum dimer models. 
A VBS order picks up a microscopic orientation that translates into a complex- valued 
order parameter in the continuum limit. Hence, we shall distinguish between the real 
(Re VBS) and imaginary (ImVBS) parts of the VBS. Triplet superconductivity (TSC) 
is also possible on bonds connecting the two sublattices. The terminology TSC will 
then also be used. To distinguish TSC with or without TRS we shall reserve the 
prefixes Re and Im for real and imaginary parts. This is a different convention for the 
use of the prefixes Re and Im than for a VBS. 

Any mass matrix that does not satisfy the SLS defined in Eq. (|22j) corresponds 
to a microscopic order parameter for which the fermion bilinear has the two lattice 
fermions sitting on the same sublattice. Microscopic examples are charge-density 
waves (CDW), spin-density waves (SDW) such as Neel ordering, orbital currents 
leading to the quantum Hall effect (QHE), spin-orbit couplings leading to the quantum 
spin Hall effect (QSHE), singlet superconductivity (SSC), or triplet superconductivity 
(TSC). 

When SU(2) spin symmetry is broken by the order parameter, we add a subindex 
x, y, or z that specifies the relevant quantization axis to the name of the mass matrix. 
Moreover, TSC with SLS must be distinguished by the 2 possible bond orientations 
(the underlying two-dimensional lattice has 2 independent vectors connecting nearest- 
neighbor sites). These 2 orientations are specified by the Pauli matrices used in the 
valley and sublattice subspaces, i.e., by the 2 pairs of numbers 02 and 32. Symmetry 
properties of all 36 PHS masses are summarized in Table [3J 

3. Majorana bound states in superconducting graphene 

The unconventional relativistic-likc band structure of graphene leads to striking 
physical phenomena, for example when graphene is placed in proximity with a 
superconductor. The problem of two dimensional Dirac fermions coupled to an s- 
wave superconductor was considered by Jackiw and Rossi [4], who have shown that 
the fermionic spectrum displays a single zero energy mode if the superconducting order 
parameter winds once about a given point in space (the center of the vortex). This 
result is to be contrasted with the case of non-relativistic s-wave superconductors, for 
which no zero mode exists in the vortex core [14] , 
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^3010 


RcVBS 


True 


True 


True 


^0020 


ImVBS 


True 


True 


True 


^3033 


CDW 


True 


True 


False 


^3003 


QHE 


False 


True 


False 


^3110 


ReVBS^ 


False 


False 


True 


x 

^0210 


ReVBS 


Th^o 1 cr> 


Fcils c 




^3310 


ReVBS 2 


False 


False 


True 


^0120 


ImVBS^ 


False 


False 


True 


x 

^3220 


TmVRCi 

J.111 V i-JkJy 


IhVi 1 CP 






^0320 


ImVBS^ 


False 


False 


True 


V 

^3103 




True 


False 


False 


^0203 


QSHE y 


True 


False 


False 


"^3303 


QSRE Z 


True 


False 


False 


^3133 


Neel x 


False 


False 


False 


^0233 


Neel y 


False 


False 


False 


^3333 


Neel z 


False 


False 


False 


^2211 


ReSSC 


True 


True 


False 


^1211 


ImSSC 


False 


True 


False 


^1002 


ReTSC 02y 


True 


False 


True 


^2002 


ImTSC 02y 


False 


False 


True 


^1102 


ReTSC 02z 


False 


False 


True 


V 

^2102 




nue 


r aise 


True 


^1302 


ReTSC 02:c 


False 


False 


True 


^2302 


ImTSC 02a , 


True 


False 


True 


^1032 


ReTSC 32y 


False 


False 


True 


^2032 


ImTSC 32!/ 


True 


False 


True 


^1132 


ReTSC 322 


True 


False 


True 


^2132 


ImTSC 32z 


False 


False 


True 


^1332 


ReTSC 32x 


True 


False 


True 


^2332 


ImTSC 32;!; 


False 


False 


True 


^1021 


ReTSC y 


True 


False 


False 


^2021 


ImTSC^ 


False 


False 


False 


^1121 


ReTSC z 


False 


False 


False 


^2121 


ImTSC z 


True 


False 


False 


^1321 


ReTSC^ 


False 


False 


False 


^2321 


ImTSC^ 


True 


False 


False 



Table 3. The 36 mass matrices with particle-hole symmetry (PHS), see Eq. 1201 . 
for the massless Dirac Hamiltonian K, from Eq. d 1 9 1) are of the form 118H and 
anticommute with /C Q . Each mass matrix can be assigned an order parameter for 
graphene. The latin subindex of the order parameter's name corresponds to the 
preferred quantization axis in SU(2) spin space. The pair of numeral subindices 
02 and 32 are used to distinguish the two unit vectors spanning two-dimensional 
space. Each mass matrix preserves or breaks time- reversal symmetry (TRS), see 
Eq. i ll' I [ ). spin-rotation symmetry (SRS), see Eq. d 2 3 D . and sublattice symmetry 
(SLS), see Eq. @2J. 



Masses and Majorana fermions in graphene 



10 



Because superconductivity mixes particles and holes, the second quantized 
operator T associated with the zero energy mode turns out to be self-adjoint, that is to 
say, r = T* . It is in the sense of being a "real" fermion that a zero mode represents a 
condensed matter realization of Majorana fermions [151 116) . An enormous amount of 
interest rests upon Majorana states regarding their possible relevance to constructing 
topological qubits: with two spatially separated vortices, each of which hosting one 
zero mode, it is possible to form a two level system (qubit) that stores the information 
non-locally. The parity of the number of zero modes per vortex core turns out to 
be fundamental in determining the stability of the qubit. If, for example, two zero 
modes exist at each vortex, generic perturbations can split those modes apart causing 
the breaking down of the stability. The general statement is then that an odd (even) 
number of zero modes per vortex implies that one can form qubits that are stable 
(unstable) against decoherence. 

Let us suppose now that graphene is brought near to a good conventional s-wave 
superconductor. By proximity effect, electron-hole pairs can tunnel between graphene 
and the superconductor. In this way, superconductivity can be induced in graphene. 
In the presence of a superconducting vortex, the low energy theory is an extension 
of that considered by Jackiw and Rossi. Zero mode solutions exist [6] but now there 
are four Majorana modes per vortex, as opposed to just an isolated Majorana mode 
in the Jackiw-Rossi system. The even number of zero modes is a direct consequence 
of the fermion doubling problem: in any TRS and local lattice theory, the number of 
Dirac cones is even [13]. Graphene can thus be regarded, as far as the discussion of 
zero modes is concerned, as four copies of the Jackiw-Rossi model. 

For the purpose, at least as a question of principle, of building topological qubits in 
graphene, or for that matter in any lattice system with Dirac particles, it is necessary 
to design a mechanism by which one can control the parity of the number of zero 
modes bound to a singly-quantized vortex in order to overcome the serious challenge 
originating from the fermion doubling problem. 

The same dilemma, but in a different context, is present in the implementation of 
lattice gauge theories, where the lattice regularization introduces spurious fermionic 
degrees of freedom. Wilson has proposed to overcome this problem by introducing 
terms "by hand" in the Hamiltonian with the effect of adding mass terms to those 
unwanted fermions, thus removing them from the low energy sector j!7| . Although 
such Wilson masses indeed remove the extra fermionic particles at tree-level, these 
perturbations have to be treated with great caution when quantum fluctuations are 
taken into account. 

We advocate that for some lattice systems considered in condensed matter physics, 
the Wilson proposal is the way to control the parity of the number of zero modes. 
Hereafter, we explain how the Wilson mechanism works for graphene. 

We recall that in graphene, electrons with spin s =j, I hop on a honeycomb lattice 
that is made of two triangular sublattices A and B. The conduction and valence bands 
touch at the two non-equivalent points K ± located at the opposite corners in the 
hexagonal first Brillouin zone (see Ref. [2] for a review). Finally, to account for the 
possibility of a SC instability, Nambu doublets are introduced with the index p and h 
to distinguish particles from their charge conjugate (holes). Hence, after linearization 
of the spectrum about the Fermi points K ± , this leads to a single-particle kinetic 
energy represented by a 16x16 dimensional matrix /C = p x a 1 + p y ct 2 - Here, a x and 
a 2 are two 16x16 dimensional Dirac matrices. 

As discussed in Sec. [2] above (and in Ref. [8]), there exists 36 distinct order 
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parameters (listed in Table such that any one, when added to /C , opens a spectral 
gap. These order parameters were identified by seeking all 16x16 matrices from the 
Clifford algebra that anticommute with /C . Among these order parameters, two (a 
real and an imaginary part) correspond to one complex-valued order parameter that 
is associated with singlet superconductivity. We shall denote the two corresponding 
16x16 matrices from the Clifford algebra by A/ RcSSC and M ImSSC and define the 
perturbation H A = A 1 Af Regsc + A 2 M ImSSC that opens the spectral gap 2|A| with 
the complex- valued A = A 1 + i A 2 parametrized by the real- valued A x and A 2 when 
added to /C . 

Next, we would like to select other order parameters among the remaining 
34(=36 — 2) masses that compete with superconductivity, i.e., masses that do not 
add in quadrature with the two superconducting masses -M Ro gsc anc ^ -^imSSC- 
Matrices corresponding to masses that add in quadrature anticommute, while matrices 
corresponding to masses that compete commute. We thus seek all 16 x 16 matrices from 
the Clifford algebra that (i) anticommute with /C and (ii) commute with 7i A . In this 
way, we find all 10 TRS-breaking order parameters listed in Table [4] that alone would 
open a gap in the Dirac spectrum if not for their competition with the gap induced 
by singlet superconductivity. Within this set of 10 matrices one can form groups of 
at most 4 matrices that are mutually commuting and therefore can be simultaneously 
diagonalized. Let us choose the 4-tuplet {ReVBS^., ImVBS y , Neel 2 , IQHE} for 
concreteness, but the results hereafter apply to any other such 4-tuplet of commuting 
mass matrices among the set of 10. Observe that any member of this 4-tuplet breaks 
TRS. This property will allow us to overcome the fermion doubling barrier [13]. It is 
then possible to show that 



U = p ■ a + A 1 M KcSSC + A 2 A/ ImSSC 

+ m i ^ReVBS x + m 2 M ImVBSa + m 3 M Ncch 

after a unitary transformation, can be brought into the form 

( n x o 
n 2 


\ 

with the 4x4 Hermitean blocks 



t/M, 



IQHE' 



H = 






^3 





\ 




H 4 ) 



(24) 



(25) 



V 



p_ 





p 

'h 
o 



j 



-p 







\ 



(26) 



whereby the complex notation p = p 1 + \p 2 is used for the momenta and x denotes 
the complex conjugate of x. It is found that 



Si 



2,3,4 



A, 



r/i = —m x + m 2 + ?7 
r/2 = na l — m 2 + m 3 
773 = m 1 + m 2 — m 3 
7/4 = — 77i j — m 2 — 77 



(27) 
(28) 
(29) 
(30) 
(31) 
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X 



3003 



IQHE 



False True False 



-^31 



3110 
0210 



3310 



ReVBS^ 
ReVBS^ 
RcVBs" 



False False True 
False False True 
False False True 



^0120 
^3220 
^0320 



ImVBS^ 
ImVBS^ 
ImVBS, 



False False True 
False False True 
False False True 



^31 

X 



3133 
0233 



NeeL, 
Neel a 
Ned', 



False False 
False False 
False False 



False 1 
False 
False 



Table 4. The 10 mass matrices with particle-hole symmetry (PHS) that 
anticommutc with a 1 and ot 2 an d commute with the singlet SC masses A^ c gsc 
and ^j m ggc Each mass matrix can be assigned an order parameter for the 
underlying microscopic model. The latin subindex of the order parameter's name 
corresponds to the preferred quantization axis in SU(2) spin space. Each mass 
matrix either preserves or breaks time-reversal symmetry (TRS), spin-rotation 
symmetry (SRS), and sublattice symmetry (SLS). An explicit representation 
defined in Ref. [8] is given in the last column. 



The breaking of the 16-dimensional matrix into four independent sectors is key to 
the ability of controlling the number of zero modes. The argument goes as follows. If 
all the rjj's are zero and the SC order parameter has a single vortex with unit winding 
number, there are 4 Majorana fermions bound to it [5J. Indeed, in this limit one 
has precisely four copies of the Jackiw-Rossi model, with any one copy delivering one 
zero-mode. 

However, as the magnitudes of the |t7 - 1 's increase, there will be a phase transition 
every time that \rjj\ = |A(r = oo)|, where |A(r = oo)| is the bulk value of the order 
parameter far away from the center of the vortex. There is no zero-mode attached 
to vortices in the j-ih copy when |t7-| > |A(r = oo)|, as can be explicitly checked. 
Indeed, this gapped phase is adiabatically connected to the limit |A(r = oo)|/|i7,-| = 
with no superconductivity, i.e., no support for a zero mode (in this gapped phase, the 
spectral symmetry of the BdG Hamiltonian prevents any change of the parity in the 
number of zero modes). Therefore, one can knock out the Majorana fermions one by 
one by taking the values of the four l^l's across the phase transitions. 

In summary, we have identified a mechanism to overcome the fermion-doubling 
barrier that can prevent the attachment of an odd number of Majorana fermions to 
the core of SC vortices in graphene-like tight-binding models. This mechanism relies 
on a topological charge that measures the parity in the number of Majorana fermions 
attached to an isolated vortex and the use of TRS-breaking order parameters that 
compete with each other and with the SC order parameter to knock out one by one 
the Majorana fermions. Therefore one can selectively switch between even and odd 
numbers. 
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